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By using different continuation methods, we unveil a wide region in the parameter space of the 
discrete cubic-quint ic complex Ginzburg-Landau equation, where several families of stable vortex 
solitons coexist. All these stationary solutions have a symmetric amplitude profile and two different 
topological charges. We also discover the dynamical formation of a variety of 'bound-state' solutions, 
composed of two or more of these vortex solitons. All of these stable composite structures persist 
in the conservative cubic limit, for high values of their power content. 



I. INTRODUCTION 

Optical beams whose phase circulates around a sin- 
gular point -or central core-, changing by 27tS times in 
each closed loop around it (with S being an integer num- 
ber), are called optical vortices. The integer number S 
is known as the topological charge of the vortex, and 
its sign defines the direction of the phase circulation. 
Usually an optical vortex has a doughnut-like shape and 
diffracts when it propagates in free space. In quantum 
information they have an enormous potential for codify- 
ing information beyond two levels using their topological 
charge value [1 . In other fields, such as biophotonics 
for example, they are useful due to their ability to affect 
the motion of particles (microorganisms) through angular 
momentum transfer [2]. Other scientific and technologi- 
cal applications for optical vortices are found in optical 
systems communication, spintronics and optical tweez- 
ers [3-[5] . These potential applications of optical vortices 
have sparkled the interest of the scientific community on 
their basic properties and characteristics. 

Diffraction is a fundamental phenomenon which leads 
to beam broadening upon propagation. In a nonlinear 
medium, self-focusing reduces diffraction whereas self- 
defocusing enhances the beam spreading. In a situation 
where the nonlinear self-focusing effects exactly balances 
diffraction, the beam can propagate as an optical spatial 
soliton, i.e. a self-trapped optical beam which preserves 
its shape upon propagation. Recently, spatial optical 
solitons have become attractive for several technologi- 
cal applications. They can be defined as self- localized 
solutions of nonlinear wave equations found in various 
physical systems [6 . Typical equations of this type in 
optics are the nonlinear Schrodinger equation (NLSE) for 
conservative systems, and the complex Ginzburg-Landau 
equation (CGLE) for dissipative ones [71|8]. The CGLE 
is a master model in which dissipative solitons [9^ are 
probably its most interesting solutions. In conserva- 
tive models, such as the ones described by the NLSE 
or its several variants, exchange of energy with the sur- 



roundings is not allowed. Self-localized solutions for the 
nonlinear Schrodinger equation originate from a balance 
between nonlinearity (e.g., the Kerr effect) and disper- 
sion/diffraction. In contrast, for dissipative systems the 
solutions also exchange energy with an external source, 
making the problem more complex and rich. In this case, 
an extra balance between gain and loss is required in or- 
der to obtain stationary solutions. In particular, dissipa- 
tive vortex solitons in continuous media have been found 
to exist for several values of and they are stable in wide 
regions of the parameter space of the CGLE [201 121]. 

In this paper we concentrate on dissipative systems 
governed by the Ginzburg-Landau equation. This model 
appears in many branches of science such as nonlinear 
optics, Bose-Einstein condensates, chemical reactions, 
super-conductivity and many others [101 HI]- Nonlin- 
ear periodic structures offer alternative ways to control 
light propagation by modifying its diffraction properties 
through the modulation of its refractive index. For in- 
stance, photonic crystals are structures of alternating re- 
fractive index that provide unprecedented control over 
light fields propagating through them; recent works show 
that lasers with square-lattices photonic crystal cavi- 
ties posses enhanced functionality and performance when 
compared to conventional lasers \TT . These systems can 
be analyzed in the framework of a set of coupled, linear 
equations which, in solid-state physics is known as the 
tight-binding approximation, while in an optics context, 
it is known as the coupled-mode approach. Stationary 
solutions obtained in this framework are called discrete 
solitons. In particular, discrete vortex solitons in con- 
servative systems have been reported on several theo- 
retical and experimental works [T6HT9] , while dissipative 
discrete solitons have been found, analytically and nu- 
merically, in one dimensional waveguide arrays [T3HT5] . 

In this paper we report the finding of a wide region 
in the parameter space of the discrete cubic-quintic com- 
plex Ginzburg-Landau equation where different discrete 
vortex solitons coexist. All the individual solutions we 
examine in this paper posses simultaneously two topo- 
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logical charges, as those reported in some of our recent 
works [22, 23^. We have studied their interactions and as 
a result the formation of bound states. 

The rest of the paper is organized as follows. In Sec- 
tion [IT] we introduce the complex cubic-quintic Ginzburg- 
Landau model that we will use in this work. Section Hill 
describes the new families of solutions we have found, 



and in Section IV we show the composite structures ob- 
tained when we let them interact as they evolve. Section 
[V| analyzes the discrete nonlinear Schrodinger equation 
case for all the solutions previously mentioned. Finally, 
section [Vll summarizes our main results and conclusions. 



II. MODEL 

Optical beam propagation in nonlinear, periodical two- 
dimensional waveguide array can be modeled by the fol- 
lowing equation: 



(1) 



Equation ([l]) is the discrete version of the complex cubic- 
quintic Ginzburg-Landau (CQGL) equation. Here, V^rn,n 
is the complex field amplitude at the (m, n) lattice site 
and iprn,n denotes its first derivative with respect to the 
propagation coordinate z. The set 

{m = -M, M}x{n = -N, N}, 

defines the array, with 2M + 1 and 2N + 1 being the 
number of sites in the horizontal and vertical directions, 
respectively. The tight binding approximation establishes 
that the field propagating in each waveguide interacts 
linearly only with nearest-neighbor fields through their 
evanescent tails. This interaction is described by the dis- 
crete diffraction operator 



m+l,n 



m,n+l 



l), 



where C is a complex parameter. Its real part denotes 
the strength of the coupling between adjacent sites and 
its imaginary part denotes the gain or loss originated by 
this coupling. The nonlinear higher-order Kerr term is 
represented by v while £ > and /i < are the co- 
efficients for cubic gain and quintic losses, respectively. 
Linear losses are accounted by a negative value of 5. 

In contrast to the conservative discrete nonlinear 
Schrodinger (DNLS) equation, the optical power, defined 
as 
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is not a conserved quantity in the present model. Nev- 
ertheless, for a self-localized solution, the power and its 
evolution will be the main quantity that we will mon- 
itor in order to identify different families of stationary 
solutions. 



We look for stationary solutions of Eq. ([T]) of the form 
'0m,n(^) = 0m,nexp(iAz) where A is real and ^m,n are 
complex amplitudes. We are interested in solutions lo- 
calized in space whose phase changes azimuthally by an 
integer number (S) of 27r along a discrete closed-circuit. 
In such a case, the self-localized solution is called a dis- 
crete vortex soliton [24 with vorticity S. By inserting 
the previous ansatz into Eq. ([T]) we obtain the following 
set of (2M + 1) X (2A^ + 1) algebraic coupled complex 
equations: 
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We solve Eq.([3| by using a multi-dimensional Newton- 
Raphson iterative algorithm. The method requires an 
initial guess, and it converges rapidly when using a 
highly-localized profile [more details can be found in 
Ref. [22]]. 



A. Linear stability analysis 

Small perturbations around the stationary solution can 
grow exponentially, leading to the destruction of the vor- 
tex soliton. A stability analysis provides us with the 
means for establishing which solutions are linearly stable. 
Let us introduce a small perturbation 0, to the localized 
stationary solution 

V^m,n = [(/>m,n + 4>mA^)]e'^^ , 4>m,n ^ C, (4) 

then, after replacing Eq.Q into Eq.Q and then lineariz- 
ing with respect to (j)^ we obtain: 

[2(1 - e)\4>^^n? + i{v - IJi)\(t>m,n\'^ " >^]4>m,n + 

[(1 - £)4L,„ + 2(J^ - ^i)\(i)m,n?<it,„]rm,n = (5) 

The solutions of the above homogeneous linear system 
can be written as 

^m,n(^) = C^,n ^Xp [7m,n^] + C^,n ^Xp [7m,n^]. (6) 

where C^'^ are integration constants and 7m, n is the dis- 
crete spectrum of the eigensystem associated with ([5|. 
The solutions are unstable if at least one eigenvalue has 
a positive real part, that is, if max{Re(7m,n)} > 0- 



III. MULTIPLICITY OF STABLE VORTEX 
SOLITON FAMILIES 

As stated above, the nonlinear gain in the system is 
mainly controlled by e; this parameter will be the only 
one that we will change in our simulations. Once we find 
a stationary solution - for a fixed set of parameters -, we 
compute its linear stability, and then change the parame- 
ters slightly and find the new solution using the previous 
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configuration. From its phase profile we can identify a 
topological charge 6* = 1 in the core - the most inner 
discrete contour - and a charge S = —3 away from the 
center. The phase profiles for families B, C, D and E 
show a topological charge 5* = — 3 in the core of these 
solutions. From B to D, the topological charge has the 
same value in the core and away from there, but the phase 
profile outside looks rotated respect to the center. For 
the last family, E, the topological charge has increased 
up to 5* = — 7 away from the center. 



IV. COMPOSITE STRUCTURES 



Figure 1 . Q versus s diagram for several vortex interconnected 
soliton families. Solid lines correspond to stable families while 
dashed lines to unstable ones. (CQGL equation parameters: 
C = 0.8, 6 = -0.9, fi = -0.1, ly = 0.1). 



one as an ansatz. In this manner we obtain all the solu- 
tion families displayed in the Q vs £ diagram shown in 
Fig.([T). The first of them (A-family), was already re- 
ported in our recent works [22l [23] . It was obtained by 
starting from the fundamental four-peaks discrete vor- 
tex soliton with = 1, after passing throughout several 
saddle-node points. A striking property of all the solu- 
tions shown in Fig.Q is that they have - simultaneously 
- two topological charges; i.e., they are two-charge vor- 
tices. Representative amplitude and phase profiles for 
these families are shown in Fig.([2|. 




Figure 2. (Color online) Color map plots for the amplitude 
(top) and phase (bottom) profiles of stable vortex solutions 
for the families marked in Fig.([T]). 

Stable families B, C, D and E - shown in Figs.Q and 
([2| - were unveiled after observing the dynamical evolu- 
tion of solutions belonging to unstable branches (dashed 
gray lines). In some cases, the stable and unstable fam- 
ilies are so close that they are nearly indistinguishable 
at the scale of Fig.Q. From the amplitude profiles we 
can see that there is a difference of four excited sites be- 
tween one stable family and the next one. Family A has 
eight main excited sites and family E has twenty four 
main peaks. All these families show a coexistence of two 
topological charges. 

The amplitude profile for case (A) shows a swirl spatial 



Next, we study the formation of bound states com- 
posed of two vortex solitons belonging to the family E 
with e = 0.9. We have chosen this family due to his 
high value of vorticity and squared symmetry equivalent 
to that of the optical lattice. We study dynamically the 
evolution of an array of two of these solutions, horizon- 
tally shifted by a certain small distance. We tested two 
initial configurations differing in their initial separation 
and, for each one of them we try a broad range of initial 
phase differences, following a procedure similar to the 
one implemented in Ref. [25 . For this purpose, we mul- 
tiply the second solution by a phase factor e*^", where 
Oa = with a = 1, 2, ...,40. A bound state is reached 
when the relative phase (A6>, defined as the phase differ- 
ence between two given sites in each solution) becomes a 
constant. In continuous and homogeneous systems, the 
separation distance also changes along the evolution and 
it becomes a constant when the bound state is formed. 
Here, in the dissipative discrete case we do not observe 
any soliton mobility and, therefore, the separation dis- 
tance remains invariant. 

We have measured the phase difference, in both con- 
figurations, for the sites enclosed by the white circles 
showed in Figjsj^a). FigsQa) and (b) show AO and Q 
versus respectively, for the first configuration shown 
in Figjsj^a). We clearly identify two attraction basins, la- 
beled as bl^ and bl^. Both (bl^) correspond to the 
lower power value shown in FigQb). This implies that 
both basins are symmetrically equivalent solutions. The 
unstable configuration is labeled as b2 and it corresponds 
to the upper power value in FigQb) [Figsjsj^c) and (d) 
show the amplitude and phase profiles of this unstable 
solution] . 

All these solutions preserve the central core structure, 
keeping the same topological charge as the initial input 
condition. Figsjsja) and (c) show the amplitude profile 
for both of them; although they are very similar, the first 
one has an extra central core (marked by a red circle), 
located at the center of the structure. For the second 
structure we can note that the column in the middle (m = 
0) is filled by small tails, without a central core. By 
taking a closed look at the rectangular contour sketched 
in Figjsj^b), we find that the phase varies continuously. 
The charge increases in the direction indicated by the 
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Figure 3. (Color online) Color map plots showing the amph- 
tude (left) and phase (right) profiles for the stable solutions 
corresponding with the basins of attraction shown in Fig[4] 
For basins bl^ the stable vortex soliton is similar to the pro- 
files shown in (a) and (b). Profiles for the b2 basin looks 
slightly different and are shown in (c) and (d). 
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Figure 4. (Color online) (a) Dynamic evolution of the relative 
phase between the two sites enclosed by the white circles in 
the vortex solutions shown in Fig.(|3|. (b) Optical power evo- 
lution for the same vortex solutions. The inset in (b) shows 
a magnification of the initial stage of evolution. 



arrows in this contour, with an accumulated charge of 
S' = 11. On the other hand, if we look how the phase 
changes along the rectangle sketched in Figjsjd), we see 
that the topological charge is not well defined on this 
contour. Indeed, the topological charge is truncated (see 
gray filled circles) meaning that this structure is not a 
composed vortex beam. Nevertheless, this profile can 
be thought as two non-interacting vortex solitons with a 
TT radians rotation between them. As we can see from 
Figjlja), any small variation in the phase leads to this 
bound solution to evolve towards the basin of attraction 
bl^. No other initial condition goes to b2, meaning 
that this is not properly a basin of attraction. So, we 
can say that vorticity is a necessary condition, achieved 
during the propagation, for the stability of this kind of 
structures. 
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Figure 5. (Color online) sm(Oi) versus / for contours (a) Fi 
and (b) Fs 



For the sake of clarity, we plot sin(^^) vs /, where / 
corresponds to the site on the inner (Fi) and outer (F2) 
discrete contours sketched in Figjsj^a). Figjsja) shows 
a good correspondence between the data (green points) 
and a sinusoidal function (gray line) with seven periods 
{S = 7) along the twenty one sites of the Fi contour. 
For the F2 contour, which contains twenty nine sites, 
we count eleven periods {S = 11) as shown in Figjsjb). 
FigjS] explicitly shows the different topological charges 
contained simultaneously in this solution. This supports 
the right identification of discrete vortex solitons, which 
is not an easy task for discrete systems. 
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Figure 6. (Color online) (a) Dynamic evolution of the rela- 
tive phase between two sites enclosed by the white circles, on 
the bound states showed in Fig.([7|. (b) optical power evo- 
lution for the same bound state. The inset in (b) shows a 
magnification of the initial stage of evolution. 

We consider now a second initial configuration where 
the two initially independent E-family vortices are placed 
closer to each other. We find a similar evolution than be- 
fore but now, there are three different stable attraction 
basins for the relative phase evolution [see Figjoja)]. Two 
of them, the lowest (b2^) and the highest (bl^), cor- 
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respond to the larger Q- value basin [b2^ in Figj6][b)]. 
The amplitude and phase profiles for these two vortices 
solutions are shown in Figsj7|^a) and (b), and (c) and 
(d), respectively. We can see from the amplitude profiles 
that these solutions lost one of the two original central 
cores (both solutions are equivalent if we perform a inver- 
sion symmetry through the n-axis). The global vorticity 
is lost, and we can see from FigslTFb) and (d) how the 
phase circulation is truncated when we move to the region 
without a core. Here, we claim that this mixed bound 
state is composed of an E-family vortex soliton and a 
staggered bright soliton (with a 7r-phase shift between 
nearest neighbors). 
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Figure 7. (Color online) Color map plots showing the ampli- 
tude (left) and phase (right) profiles for the stable solutions 
corresponding with the basins of attraction showed in Fig[6] 
For basin b2^ the stable structure is similar to the profiles 
shown in (a) and (b), while (c) and (d) correspond with b2^. 
For the remaining hi basin, the vortex soliton is similar to the 
profiles shown in (e) and (f). 

The third basin (bl), which corresponds to the lower 
Q-value basin in Figj6[b), has the amplitude and phase 
profiles displayed in Figs|7|^e) and (f). We clearly see that 
it preserves the initial two central cores keeping the same 
topological charge as the initial condition. Unlike the 
previous two basins, the global topological charge of this 
solution is well-defined. As for the first configuration, 
there are also two different topological charges for this 
composite vortex soliton. Again, to corroborate this, we 



plot sin(6>/) vs / in order to show in detail the topological 
charge along two different contours. The first one. Fa, 
corresponds to the inner rectangular contour sketched in 
Figj7|^e), while F4 corresponds with the outer rectangular 
contour sketched in the same figure. Figjsja) shows how 
the inner charge is 5 = 7 while Figjsj^b) indicates a charge 
= 11 for contour F4. 




1 4 7 10 13 16 19 22 25 
/ 

Figure 8. (Color online) shi{Oi) vs / (site number) diagram 
for the first (a) and second (b) discrete contours shown on 
the Fig|6|a) for the bound state vortex soliton. 

We show now two more examples of composite struc- 
tures built from an initial superposition of two-and four 
solutions taken from the D and E families shown in Fig{l] 
In both cases, the values of the parameters used are: 
C = 0.8, 6 = -0.9, e = 0.9, /i = -0.1 and u = 0.1. 
The typical propagation distance was z ^ 300, enough 
for the power content to become constant. 

Figure [9] shows three stable solutions obtained by su- 
perposing two vortex solitons belonging to the D-family 
m Figjl] The first one is constructed by overlapping two 
of these vortices spaced by one site between their central 
cores. Figsj9|a) and (b) show the amplitude and phase 
profiles for this stable solution. We note that this state 
has only one central core, located halfway between the 
initial ones. On the other hand, the phase profile shows 
a charge 6* = 5 at the inner contour and rotated with 
respect to the next discrete contours. 

The next configuration is constructed in the same man- 
ner as the previous one but now, the center-to-center dis- 
tance between the cores has been increased to two sites. 
Figure [9]^c) shows a dynamically stable solution with two 
central cores located at the same positions as the ini- 
tial condition. The phase profile [see Figj9]^d)] shows a 
value of 5 = 5 for the topological charge, as in the previ- 
ous case. A third stable composed structure is obtained 
by superposing again two vortex solitons with an initial 
center-to-center distance of three sites. The amplitude 
profile for the new dynamically stable solution has one 
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Figure 9. (Color online) Color map plots showing the ampli- 
tude (left) and phase (right) profiles for three stable vortex 
solutions obtained by superposing two vortex solitons belong- 
ing to the D family, at different initial distances. 



horizontally elongated core, along two lattice sites, as 
shown in Figj9|e), and its topological charge has two dif- 
ferent values as shown in Figjoj^f). Indeed, the innermost 
discrete contour exhibits a charge 5 = 6, while the re- 
maining contours have a charge 6* = 10. 



Finally, we show another example of a composed struc- 
ture that was obtained by combining four solutions be- 
longing to the E family . We locate each E-vortex by 
placing their central cores forming the vertices of a 8 x 8 
square. We use this configuration as initial condition for 
model ([1]) and find a dynamically stable stationary solu- 
tion. Figs. [lOj^a) and (b) show the amplitude and phase 
profiles for this composite solution. We observe a spatial 
amplitude distribution similar to the initial condition, 
where the four initial cores preserve their initial position 
and vorticity. In addition, an extra phase core appears 
at the lattice center (n, m = 0), around which a 5* = — 1 
topological charge can be observed. If we follow a new 
contour that encloses all the sites with large amplitude, 
the topological charge measured will be S = 15. 
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Figure 10. (Color online) Color map plots for the amplitude 
(a) and phase (b) profiles for a dynamically stable solution ob- 
tained by combining four solutions belonging to the E family 



V. SCHRODINGER LIMIT 

Most of the present experiments with optical beams are 
performed under conditions that closely meet the cubic 
conservative limit. So, we are interested in knowing if 
all these previous dissipative structures can be observed 
also here. In this scenario, all the dissipative parameters 
are suppressed, i.e., u = /i = s = S = 0, and model ([3| 
reduces to the discrete NLSE equation 



0. 



(7) 



Taking as an ansatz for the decoupled limit for some of 
the solutions previously described, we have obtained the 
same kind of structures for the discrete NLSE. We have 
constructed the corresponding families and also com- 
puted their corresponding stability region. We note that 
all these solutions exist in the Schrodinger limit, but are 
stable only for high values of their optical power content. 

In Fig. ( 11 ) we show five families of stationary solutions 
of the NLSE that correspond with some of the solutions 
reported in the previous sections. Namely, the blue line 
is the family corresponding with the solution displayed 
in Fig.([3|i-b), the black line with the Fig.([7^-d), and the 
red, green and gray line with Fig.([9^-b), Fig.([9]3-d) and 
Fig.([9^-f) respectively. The stable (unstable) solutions 
are represented by continuous (dashed) lines. The inset 



located in the upper left corner at Fig.([ll| shows a mag- 
nified view of the region close to the 
zone). As usual [22 



inear band (gray 
each one of these families tends 
to increase its power steeply after passing through the 
point of minimum optical power. The other inset at the 
lower right corner shows a zoom of the black and blue 
curves, which are almost indistinguishable because they 
have very similar spatial profiles. 



VI. SUMMARY AND CONCLUSIONS 

We have reported several new families of discrete vor- 
tex solitons, characterized for having two topological 
charges simultaneously, and coexisting for the same set 
of parameters. By superposing two or more of these vor- 
tices, we have been able to produce new, dynamically 
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Figure 11. Q versus A diagram displaying the families cor- 
responding with some of dissipative solutions reported previ- 
ously, but for the conservative cubic case. 



stable composite vortex solitons that are also endowed 
with multiple vorticity charges. Additionally, we have 
shown that these composite structures persist in the con- 
servative limit and they are stable for high values of their 
power content. 
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